%code for calibration of parameters and calculation of mcf by using the model
%with informality

clear all
clc
syms thetai gamma thetaf tauk taun tauc alpha fgdp igdp phi nrate b delta

filename = 'data.xlsx';
data = xlsread(filename);
rho=0; %enforcement parameter is 0 
Time=16*330;
taui=0; %there is no tax for informal sector 
[r,c]=size(data);

 for i=1:r
    
    alpha=data(i,1);
    b=data(i,2);
    delta=data(i,3);
    
    tauk=data(i,4);
    taun=data(i,5);
    tauc=data(i,6);
    
    fgdp=data(i,7);
    igdp=data(i,8);
    phi=data(i,9);
    nrate=data(i,10);
   
labori1=(((1-rho*taui)*gamma*thetai)/((1-taun)*(1-alpha)*thetaf))^(1/(1-gamma));
labori2=(((1/b-1)+delta)/(alpha*(1-tauk)*thetaf))^(alpha/((1-alpha)*(1-gamma)));
labori=labori1*labori2;
laborf1=(Time-labori)*gamma*(1-rho*taui)*thetai*labori^(gamma-1)-phi*(1-rho*taui)*thetai*labori^gamma;
laborf2a=gamma*(1*-rho*taui)*thetai*labori^(gamma-1);
laborf2b=phi*[(alpha*(1-tauk)+(1-alpha)*(1-taun))*thetaf*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(alpha/(1-alpha))-delta*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(1/(1-alpha))];
laborf=laborf1/(laborf2a+laborf2b);
K1=(((1-tauk)*thetaf*alpha)/(1/b-1+delta))^(1/(1-alpha));
K=laborf*K1;
Yf=thetaf*(K^(alpha))*(laborf^(1-alpha));
Yi=thetai*(labori^(gamma));
C=((1-tauk)*alpha*Yf+(1-taun)*(1-alpha)*Yf+(1-rho*taui)*Yi-delta*K)/(1+tauc);
infrate=Yi/Yf;
lrate=labori/laborf;
leisure=Time-laborf-labori;
R=tauc*C+tauk*alpha*Yf+taun*(1-alpha)*Yf+rho*taui*Yi;

%calibration of parameters -thetai, tetaf, gamma-
z0=[0.5,100,0.5];
options= optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt',...
    'MaxFunctionEvaluations',1500);

res=lsqnonlin(@(z) funcs (z,tauk,taun,b,alpha,fgdp,igdp, phi ,nrate,delta),z0, [0,0,0],[],options);

%calculation of the values by using calibrated varaibles
Rres=subs(R,{thetaf,thetai,gamma},{res(1),res(2),res(3)});
leisureres=subs(leisure,{thetaf,thetai,gamma},{res(1),res(2),res(3)});
laborires=subs(labori,{thetaf,thetai,gamma},{res(1),res(2),res(3)});
laborfres=subs(laborf,{thetaf,thetai,gamma},{res(1),res(2),res(3)});
Kres=subs(K,{thetaf,thetai,gamma},{res(1),res(2),res(3)});
Cres=subs(C,{thetaf,thetai,gamma},{res(1),res(2),res(3)});
Yires=subs(Yi,{thetaf,thetai,gamma},{res(1),res(2),res(3)});
Yfres=subs(Yf,{thetaf,thetai,gamma},{res(1),res(2),res(3)});

data(i,11)=res(1);
data(i,12)=res(2);
data(i,13)=res(3);

data (i,14)=eval(leisureres);
data (i,15)=eval(laborires);
data (i,16)=eval(laborfres);
data (i,17)=eval(Kres);
data (i,18)=eval(Cres);
data (i,19)=eval(Yfres);
data (i,20)=eval(Yires);
data (i,21)=eval(Rres);
 end
 
  
 filename='input.xlsx';
 xlswrite(filename,data);
clear all
clc


filename = 'input.xlsx';
mcfdata = xlsread(filename);
rho=0;
Time=16*330;
taui=0;

syms tauk taun tauc thetai gamma thetaf alpha b phi delta
labori1=(((1-rho*taui)*gamma*thetai)/((1-taun)*(1-alpha)*thetaf))^(1/(1-gamma));
labori2=(((1/b-1)+delta)/(alpha*(1-tauk)*thetaf))^(alpha/((1-alpha)*(1-gamma)));
labori=labori1*labori2;
laborf1=(Time-labori)*gamma*(1-rho*taui)*thetai*labori^(gamma-1)-phi*(1-rho*taui)*thetai*labori^gamma;
laborf2a=gamma*(1*-rho*taui)*thetai*labori^(gamma-1);
laborf2b=phi*[(alpha*(1-tauk)+(1-alpha)*(1-taun))*thetaf*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(alpha/(1-alpha))-delta*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(1/(1-alpha))];
laborf=laborf1/(laborf2a+laborf2b);
K1=(((1-tauk)*thetaf*alpha)/(1/b-1+delta))^(1/(1-alpha));
K=laborf*K1;
Yf=thetaf*(K^(alpha))*(laborf^(1-alpha));
Yi=thetai*(labori^(gamma));

[r,c]=size(mcfdata);

C=((1-tauk)*alpha*Yf+(1-taun)*(1-alpha)*Yf+(1-rho*taui)*Yi-delta*K)/(1+tauc);
R=tauc*C+tauk*alpha*Yf+taun*(1-alpha)*Yf+rho*taui*Yi;
leisure=Time-laborf-labori;
utility=log(C)+phi*log(leisure);
turevRtaui=diff(R,taui);
turevUtaui=diff(utility,taui);
turevRtauc=diff(R,tauc);
turevUtauc=diff(utility,tauc);
mcpftaui=-turevUtaui/turevRtaui;
turevRtauk=diff(R,tauk);
turevUtauk=diff(utility,tauk);
mcpftauk=-(turevUtauk/turevRtauk);
turevRtaun=diff(R,taun);
turevUtaun=diff(utility,taun);
mcpftaun=-(turevUtaun/turevRtaun);
mcpftauc=-(turevUtauc/turevRtauc);


for i=1:r
  %calculation of marginal cost of public funds
    mcfk=subs(mcpftauk,{alpha,b,delta,tauk,taun,tauc,phi,thetaf,thetai,gamma},{mcfdata(i,1),mcfdata(i,2),mcfdata(i,3),mcfdata(i,4),mcfdata(i,5),mcfdata(i,6),mcfdata(i,9),mcfdata(i,11),mcfdata(i,12),mcfdata(i,13)});
    mcfc=subs(mcpftauc,{alpha,b,delta,tauk,taun,tauc,phi,thetaf,thetai,gamma},{mcfdata(i,1),mcfdata(i,2),mcfdata(i,3),mcfdata(i,4),mcfdata(i,5),mcfdata(i,6),mcfdata(i,9),mcfdata(i,11),mcfdata(i,12),mcfdata(i,13)});
    mcfn=subs(mcpftaun,{alpha,b,delta,tauk,taun,tauc,phi,thetaf,thetai,gamma},{mcfdata(i,1),mcfdata(i,2),mcfdata(i,3),mcfdata(i,4),mcfdata(i,5),mcfdata(i,6),mcfdata(i,9),mcfdata(i,11),mcfdata(i,12),mcfdata(i,13)});
    
    mcfdata(i,22)=eval(mcfk);
    mcfdata(i,23)=eval(mcfc);
    mcfdata(i,24)=eval(mcfn);
 
end

filename = 'output.xlsx';
xlswrite(filename,mcfdata);
